Made maps and graphs of our models
Altered our models to get them to converge.
Predicting hotspots of specific crimes and how they changed over time and space.
Are certain types of crimes more prevalent in certain areas than others?
How have these hotspots changed over time? If they have changed over time, how has it changed relative to changes in income, demographics, age, etc?
Examining police traffic stops and subsequent arrests to determine whether certain ethnic minorities are stopped/arrested more than others.
Are certain ethnic groups more likely to be arrested/stopped in certain areas than others? For instance, are blacks/hispanics more likely to be stopped/arrested in predominantly white neighborhoods?
We would like to examine the relationship between traffic stops/arrests and the actual crime rates of different ethnic groups. For instance, are blacks/hispanics getting stopped/arrested at a disproportionately greater rate than they are committing crimes?
library(tidyr)
library(rjags)
library(ggplot2)
library(dplyr)
library(tidyverse)
library(MacBayes)
library(choroplethrZip)
library(rgdal)
cleanup <- theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_rect(fill = 'white', colour = 'white'),
axis.line = element_line(colour = "white"),
axis.ticks=element_blank(), axis.text.x=element_blank(),
axis.text.y=element_blank())
#Load data
Crime = read.csv('Crime_sample_2010_to_2017.csv')
Crime = Crime %>%
mutate(Date = as.POSIXct(strptime(Date, "%Y-%m-%d %H:%M:%S"))) %>%
filter(!is.na(Latitude))
## Warning in strptime(Date, "%Y-%m-%d %H:%M:%S"): unknown timezone 'zone/tz/
## 2017c.1.0/zoneinfo/America/Chicago'
#load shapefile
zip_shp <- readOGR(dsn="Boundaries - ZIP Codes", layer="geo_export_e1262361-5c82-45ee-8427-ef228e06dc4a")
## OGR data source with driver: ESRI Shapefile
## Source: "Boundaries - ZIP Codes", layer: "geo_export_e1262361-5c82-45ee-8427-ef228e06dc4a"
## with 61 features
## It has 4 fields
temp <- zip_shp
zip_df <- fortify(temp, region = "zip")
#needs to reassign CRS for shapefile
new_CRS <- CRS("+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
zip_shp <- sp::spTransform(zip_shp, new_CRS)
#summarize crime count per year
years <- unique(Crime$Year)
Hotspot_Crime <- data.frame()
for(i in 1:length(years)){
year <- years[i]
Current_Year <- Crime %>% filter(Year == year)
locations <- with(Current_Year, as.data.frame(cbind(Longitude, Latitude)))
coordinates(locations) <- ~Longitude+Latitude
proj4string(locations) <- CRS("+init=epsg:4326")
#prepping the data: project points to polygons
proj4string(zip_shp) <- new_CRS
proj4string(locations) <- new_CRS
by_zip <- over(locations, zip_shp)
by_zip <- by_zip %>%
group_by(zip) %>%
dplyr::summarise(total=n()) %>%
filter(!is.na(zip)) %>%
mutate(id = as.character(zip))
total_map <- left_join(zip_df, by_zip)
total_map <- total_map %>% mutate(Year = year)
Hotspot_Crime <- rbind(Hotspot_Crime, total_map)
}
data("df_zip_demographics")
Demographics <- df_zip_demographics %>% mutate(zip = region)
Hotspot_Crime <- inner_join(Hotspot_Crime, Demographics, by="zip")
## Warning: Column `zip` joining factor and character vector, coercing into
## character vector
Model_Data <- Hotspot_Crime %>%
select(-c(lat, long, order, piece, group, hole)) %>% distinct
Crime_count <- Model_Data %>%
select(Year, zip, total) %>%
spread(key=Year, value=total)
rownames(Crime_count) <- Crime_count$zip
Crime_count <- Crime_count %>% select(-zip)
Crime_count[is.na(Crime_count)] <- 0
head(Crime_count)
## 2010 2011 2012 2013 2014 2015 2016 2017
## 60601 8 6 5 8 7 6 7 11
## 60602 4 8 1 4 5 10 8 9
## 60603 1 5 3 2 1 2 9 4
## 60604 1 6 3 1 4 3 1 3
## 60605 10 7 13 6 9 10 10 8
## 60606 2 6 4 3 1 4 3 2
Population <- Model_Data %>%
select(total_population, zip, Year) %>%
spread(key=Year, value=total_population) %>%
select('2010', zip)
rownames(Population) <- Population$zip
Population <- Population %>% select(-zip)
Population[is.na(Population)] <- 0
head(Population)
## 2010
## 60601 9732
## 60602 1463
## 60603 882
## 60604 415
## 60605 24972
## 60606 2789
Population.vector <- Population$`2010`
time <- sort(unique(Model_Data$Year))
percent_nonwhite <- 100 - Model_Data$percent_white
Plot of \(log(y_{ij}/n_{i})\) vs year
Crime_count_log_calc <- log(Crime_count/Population.vector)
Crime_count_log <- Crime_count_log_calc %>%
rownames_to_column() %>%
gather(key=year, value=crime_count_log, -rowname) %>%
mutate(zip=rowname) %>%
select(zip, year, crime_count_log) %>%
mutate(year=as.numeric(year))
head(Crime_count_log)
## zip year crime_count_log
## 1 60601 2010 -7.103733
## 2 60602 2010 -5.901950
## 3 60603 2010 -6.782192
## 4 60604 2010 -6.028279
## 5 60605 2010 -7.822925
## 6 60606 2010 -7.240291
ggplot(Crime_count_log, aes(x=year, y=crime_count_log, colour=zip)) + geom_point() + geom_line()
Plot of \(log(y_{ij}/n_{i})\) vs per capita income
Model_Data_1 <- Model_Data %>%
filter(Year == 2010) %>%
mutate(percent_nonwhite = 100 - percent_white) %>%
select(zip, per_capita_income, percent_nonwhite) %>%
right_join(Crime_count_log, by="zip")
head(Model_Data_1)
## zip per_capita_income percent_nonwhite year crime_count_log
## 1 60601 89136 37 2010 -7.103733
## 2 60602 66748 38 2010 -5.901950
## 3 60603 73239 39 2010 -6.782192
## 4 60604 135807 20 2010 -6.028279
## 5 60605 58338 42 2010 -7.822925
## 6 60606 83005 25 2010 -7.240291
ggplot(Model_Data_1, aes(x=log(per_capita_income), y=crime_count_log, colour=percent_nonwhite)) + geom_point() + labs(x="Log of Per Capita Income", y="Log of Crime Rate", title="Relationship between Per Capita Income and Log of Crime Rate")
Plot of \(log(y_{ij}/n_{i})\) vs nonwhite percentage
ggplot(Model_Data_1, aes(x=log(percent_nonwhite), y=crime_count_log)) + geom_point() + labs(x="Log of Nonwhite Percentage", y="Log of Crime Rate", title="Relationship between Nonwhite Percentage and Log of Crime Rate")
\[y_{ij} \sim Bin(p_{ij}, n_{i})\] \[logit(p_{ij}) = \beta * time_{j}\]
\[\beta \sim N(0, 1000^2)\]
i: zipcodes in Chicago (58 total in our data)
j: years (from 2010 to 2017)
y[i, j]: observed count of violent crimes in Chicago by year and zipcode
p[i, j]: probability of being a victim in location i and time j
n[i]: the population in zipcode i
\(\beta\): time trend
time[j]: the time in years
#Data
Crime_count <- Model_Data %>%
group_by(Year) %>%
dplyr::summarise(total.stop = sum(total))
head(Crime_count)
## # A tibble: 6 x 2
## Year total.stop
## <int> <int>
## 1 2010 1474
## 2 2011 1448
## 3 2012 1340
## 4 2013 1258
## 5 2014 1090
## 6 2015 1117
#Specify the model
crime_model_0 <- "model{
#Data
for(j in 1:8){
y[j] ~ dbin(p[j], n)
logit(p[j]) = beta*time[j]
}
#Priors
beta ~ dnorm(0, 1/100^2)
}"
#set up an algorithm to simulate the posterior by
#combining the model (nba_model) and data (y)
#set the random number seed
crime_jags_0 <- jags.model(textConnection(crime_model_0), data=list(y=Crime_count$total.stop, n=sum(Population.vector), time=Crime_count$Year),inits=list(.RNG.name="base::Wichmann-Hill", .RNG.seed=1989))
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 8
## Unobserved stochastic nodes: 1
## Total graph size: 40
##
## Initializing model
#simulate a sample from the posterior
#note that we specify both mu and tau variables
crime_sim_0 <- coda.samples(crime_jags_0, variable.names=c("p", "beta"), n.iter=50000)
#store the samples in a data frame:
crime_data_0 <- data.frame(crime_sim_0[[1]])
running_mean_plot(crime_data_0$beta)
ggplot(crime_data_0, aes(x=beta)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Effect of Time", x="Time Effect", y="Density")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(crime_data_0, aes(x=p.1.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2010", x="Probability of Crime", y="Density")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(crime_data_0, aes(x=p.2.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2011", x="Probability of Crime", y="Density")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(crime_data_0, aes(x=p.3.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2012", x="Probability of Crime", y="Density")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(crime_data_0, aes(x=p.4.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2013", x="Probability of Crime", y="Density")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(crime_data_0, aes(x=p.5.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2014", x="Probability of Crime", y="Density")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(crime_data_0, aes(x=p.6.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2015", x="Probability of Crime", y="Density")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(crime_data_0, aes(x=p.7.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2016", x="Probability of Crime", y="Density")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(crime_data_0, aes(x=p.8.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2017", x="Probability of Crime", y="Density")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
crime_model_0_CI <- matrix(nrow=9, ncol=2)
for (i in 1:9) {
crime_model_0_CI[i,] <- quantile(crime_data_0[[i]], c(0.025, 0.975))
}
crime_model_0_mean <- matrix(nrow=9, ncol=1)
for (i in 1:9) {
crime_model_0_mean[i, ] <- mean(crime_data_0[[i]])
}
crime_model_0_table <- cbind(crime_model_0_CI, crime_model_0_mean)
rownames(crime_model_0_table) <- names(crime_data_0)
colnames(crime_model_0_table) <- c("2.5%", "97.5%", "Mean")
as.table(crime_model_0_table)
## 2.5% 97.5% Mean
## beta -0.0038479706 -0.0038283390 -0.0038380756
## p.1. 0.0004373144 0.0004549075 0.0004461179
## p.2. 0.0004356356 0.0004531701 0.0004444097
## p.3. 0.0004339632 0.0004514393 0.0004427081
## p.4. 0.0004322972 0.0004497151 0.0004410130
## p.5. 0.0004306377 0.0004479975 0.0004393244
## p.6. 0.0004289845 0.0004462865 0.0004376422
## p.7. 0.0004273377 0.0004445820 0.0004359665
## p.8. 0.0004256971 0.0004428839 0.0004342971
i: zipcodes in Chicago (58 total in our data)
j: years (from 2010 to 2017)
y[i, j]: observed count of violent crimes in Chicago by year and zipcode
p[i, j]: probability of being a victim in location i and time j
n[i, j]: the population in zipcode i in time j.
time[j]: the time in years
\(\delta\)[i]: the spatial temporal interaction term
We are modeling the count of crime incidents by year and zipcodes using a binomial distribution, in which the parameters are the probability of being a victim and the total population in location i and year j.
We modeled the probability of being a victim using a logit function of the grand crime rate of Chicago and the inclusion of several random effects. It is also a function of time and its interaction with location.
\[y_{ij} \sim Bin(p_{ij}, n_{i})\] \[logit(p_{ij}) = \delta_{i} * time_{j}\]
\[\delta_{i} \sim N(0, prec.delta)\]
\[prec.delta \sim Gamma(0.5, 0.0005)\]
#Data
Crime_count <- Model_Data %>%
select(Year, zip, total) %>%
spread(key=Year, value=total)
rownames(Crime_count) <- Crime_count$zip
Crime_count <- Crime_count %>% select(-zip)
Crime_count[is.na(Crime_count)] <- 0
head(Crime_count)
## 2010 2011 2012 2013 2014 2015 2016 2017
## 60601 8 6 5 8 7 6 7 11
## 60602 4 8 1 4 5 10 8 9
## 60603 1 5 3 2 1 2 9 4
## 60604 1 6 3 1 4 3 1 3
## 60605 10 7 13 6 9 10 10 8
## 60606 2 6 4 3 1 4 3 2
#Specify the model
crime_model_1 <- "model{
#Data
for(i in 1:58) {
for(j in 1:8){
y[i,j] ~ dbin(p[i,j], n[i])
logit(p[i,j]) = delta[i]*time[j]
}
}
#Priors
for(i in 1:58){
delta[i] ~ dnorm(0, prec.delta)
}
#Hyperpriors
prec.delta ~ dgamma(0.5, 0.0005)
}"
#set up an algorithm to simulate the posterior by
#combining the model (nba_model) and data (y)
#set the random number seed
crime_jags_1 <- jags.model(textConnection(crime_model_1), data=list(y=Crime_count, n=Population.vector, time=time),
inits=list(.RNG.name="base::Wichmann-Hill", .RNG.seed=1989))
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 464
## Unobserved stochastic nodes: 59
## Total graph size: 1520
##
## Initializing model
#simulate a sample from the posterior
#note that we specify both mu and tau variables
crime_sim_1 <- coda.samples(crime_jags_1, variable.names=c("p", "delta"), n.iter=60000)
#store the samples in a data frame:
crime_data_1 <- data.frame(crime_sim_1[[1]])
ggplot(crime_data_1, aes(x=delta.1.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Crime Differential Growth in 2010", x="Crime Differential Growth", y="Density")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(crime_data_1, aes(x=p.1.1.)) + geom_histogram(aes(y=..density..), colour="white") + labs(title="Posterior Distribution of the Probability of Crime in 2010 in Zipcode 1", x="Probability of Crime", y="Density")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
crime_model_1_CI <- matrix(nrow=length(names(crime_data_1)), ncol=2)
for (i in 1:length(names(crime_data_1))) {
crime_model_1_CI[i, ] <- quantile(crime_data_1[[i]], c(0.025, 0.975))
}
crime_model_1_mean <- matrix(nrow=length(names(crime_data_1)), ncol=1)
for (i in 1:length(names(crime_data_1))) {
crime_model_1_mean[i, ] <- mean(crime_data_1[[i]])
}
crime_model_1_table <- cbind(crime_model_1_CI, crime_model_1_mean)
rownames(crime_model_1_table) <- names(crime_data_1)
colnames(crime_model_1_table) <- c("2.5%", "97.5%", "Mean")
as.table(crime_model_1_table)
## 2.5% 97.5% Mean
## delta.1. -3.711729e-03 -3.456400e-03 -3.580096e-03
## delta.2. -2.867138e-03 -2.586290e-03 -2.722396e-03
## delta.3. -2.969974e-03 -2.586799e-03 -2.769938e-03
## delta.4. -2.718293e-03 -2.298821e-03 -2.498395e-03
## delta.5. -4.051228e-03 -3.821873e-03 -3.933484e-03
## delta.6. -3.588951e-03 -3.197713e-03 -3.382634e-03
## delta.7. -3.879030e-03 -3.684672e-03 -3.779524e-03
## delta.8. -4.119015e-03 -3.975270e-03 -4.046257e-03
## delta.9. -3.750831e-03 -3.637632e-03 -3.693651e-03
## delta.10. -4.088548e-03 -3.894629e-03 -3.989673e-03
## delta.11. -3.756364e-03 -3.595060e-03 -3.674093e-03
## delta.12. -3.611611e-03 -3.482366e-03 -3.545865e-03
## delta.13. -4.320107e-03 -4.105318e-03 -4.209816e-03
## delta.14. -4.097897e-03 -3.946837e-03 -4.021145e-03
## delta.15. -3.980550e-03 -3.808793e-03 -3.892875e-03
## delta.16. -4.056067e-03 -3.886810e-03 -3.970244e-03
## delta.17. -3.837891e-03 -3.730789e-03 -3.783815e-03
## delta.18. -4.166437e-03 -4.030778e-03 -4.097490e-03
## delta.19. -3.668010e-03 -3.564934e-03 -3.615850e-03
## delta.20. -3.587042e-03 -3.497329e-03 -3.541414e-03
## delta.21. -3.400624e-03 -3.292082e-03 -3.346049e-03
## delta.22. -3.974423e-03 -3.823658e-03 -3.897974e-03
## delta.23. -3.749702e-03 -3.652031e-03 -3.700433e-03
## delta.24. -3.351117e-03 -3.254202e-03 -3.302190e-03
## delta.25. -4.264045e-03 -4.100337e-03 -4.180989e-03
## delta.26. -4.122069e-03 -3.945685e-03 -4.032168e-03
## delta.27. -3.631162e-03 -3.536932e-03 -3.583492e-03
## delta.28. -3.966688e-03 -3.860937e-03 -3.913238e-03
## delta.29. -4.278376e-03 -4.086159e-03 -4.179781e-03
## delta.30. -4.747881e-03 -4.358148e-03 -4.543484e-03
## delta.31. -4.150986e-03 -4.012304e-03 -4.080392e-03
## delta.32. -4.613369e-03 -4.136416e-03 -4.360774e-03
## delta.33. -4.320027e-03 -4.140211e-03 -4.228669e-03
## delta.34. -3.480584e-03 -3.371023e-03 -3.425328e-03
## delta.35. -3.604060e-03 -3.491357e-03 -3.546880e-03
## delta.36. -4.438547e-03 -4.213139e-03 -4.321765e-03
## delta.37. -4.011126e-03 -3.890165e-03 -3.949733e-03
## delta.38. -4.046383e-03 -3.898156e-03 -3.970756e-03
## delta.39. -4.120865e-03 -3.969203e-03 -4.043773e-03
## delta.40. -3.825775e-03 -3.611044e-03 -3.715930e-03
## delta.41. -3.908333e-03 -3.764577e-03 -3.835530e-03
## delta.42. -3.463873e-03 -3.365678e-03 -3.414195e-03
## delta.43. -4.320507e-03 -4.100282e-03 -4.207425e-03
## delta.44. -4.819598e-03 -4.399923e-03 -4.598882e-03
## delta.45. -4.029310e-03 -3.903263e-03 -3.965421e-03
## delta.46. -3.567763e-03 -3.456856e-03 -3.511583e-03
## delta.47. -3.645227e-03 -3.541671e-03 -3.592990e-03
## delta.48. -4.309039e-03 -4.080336e-03 -4.191130e-03
## delta.49. -3.715558e-03 -3.564806e-03 -3.638507e-03
## delta.50. -3.663092e-03 -3.465155e-03 -3.561759e-03
## delta.51. -4.471091e-03 -4.157674e-03 -4.307555e-03
## delta.52. -4.593938e-03 -4.240701e-03 -4.410293e-03
## delta.53. -4.154460e-03 -3.995037e-03 -4.072562e-03
## delta.54. -4.328963e-03 -4.090126e-03 -4.206603e-03
## delta.55. -4.244785e-03 -4.029039e-03 -4.134586e-03
## delta.56. -3.888913e-03 -3.554766e-03 -3.714995e-03
## delta.57. -4.794647e-03 -4.448463e-03 -4.614012e-03
## delta.58. -4.593432e-03 -4.245760e-03 -4.413045e-03
## p.1.1. 5.749939e-04 9.602428e-04 7.554488e-04
## p.2.1. 3.131995e-03 5.494800e-03 4.227676e-03
## p.3.1. 2.548631e-03 5.489217e-03 3.877168e-03
## p.4.1. 4.219669e-03 9.750573e-03 6.699385e-03
## p.5.1. 2.906884e-04 4.608558e-04 3.708225e-04
## p.6.1. 7.358186e-04 1.614035e-03 1.136081e-03
## p.7.1. 4.108583e-04 6.071116e-04 5.042648e-04
## p.8.1. 2.536708e-04 3.386202e-04 2.943978e-04
## p.9.1. 5.315555e-04 6.672752e-04 5.972580e-04
## p.10.1. 2.696861e-04 3.981812e-04 3.305846e-04
## p.11.1. 5.256803e-04 7.268444e-04 6.222627e-04
## p.12.1. 7.030785e-04 9.114562e-04 8.040822e-04
## p.13.1. 1.693430e-04 2.607494e-04 2.126423e-04
## p.14.1. 2.646669e-04 3.585291e-04 3.097355e-04
## p.15.1. 3.350469e-04 4.731263e-04 4.011305e-04
## p.16.1. 2.878756e-04 4.044865e-04 3.433503e-04
## p.17.1. 4.462611e-04 5.533940e-04 4.982303e-04
## p.18.1. 2.306134e-04 3.028821e-04 2.655262e-04
## p.19.1. 6.277760e-04 7.721824e-04 6.980831e-04
## p.20.1. 7.386444e-04 8.844757e-04 8.103880e-04
## p.21.1. 1.074039e-03 1.335537e-03 1.200276e-03
## p.22.1. 3.391972e-04 4.592063e-04 3.966882e-04
## p.23.1. 5.327634e-04 6.482515e-04 5.889284e-04
## p.24.1. 1.186280e-03 1.441042e-03 1.310311e-03
## p.25.1. 1.895382e-04 2.633725e-04 2.247577e-04
## p.26.1. 2.521184e-04 3.593598e-04 3.032732e-04
## p.27.1. 6.760041e-04 8.168536e-04 7.448089e-04
## p.28.1. 3.445098e-04 4.260690e-04 3.841264e-04
## p.29.1. 1.841574e-04 2.709842e-04 2.256018e-04
## p.30.1. 7.167902e-05 1.568791e-04 1.102606e-04
## p.31.1. 2.378859e-04 3.143366e-04 2.748372e-04
## p.32.1. 9.392947e-05 2.449535e-04 1.606904e-04
## p.33.1. 1.693701e-04 2.430927e-04 2.043642e-04
## p.34.1. 9.147224e-04 1.139806e-03 1.023656e-03
## p.35.1. 7.138241e-04 8.951459e-04 8.020259e-04
## p.36.1. 1.334728e-04 2.099544e-04 1.699001e-04
## p.37.1. 3.150820e-04 4.017686e-04 3.571354e-04
## p.38.1. 2.935322e-04 3.953701e-04 3.427030e-04
## p.39.1. 2.527294e-04 3.427737e-04 2.959657e-04
## p.40.1. 4.572566e-04 7.038798e-04 5.736094e-04
## p.41.1. 3.873673e-04 5.170781e-04 4.495990e-04
## p.42.1. 9.459398e-04 1.152103e-03 1.046489e-03
## p.43.1. 1.692068e-04 2.634018e-04 2.137256e-04
## p.44.1. 6.205724e-05 1.442461e-04 9.893437e-05
## p.45.1. 3.037772e-04 3.913338e-04 3.461138e-04
## p.46.1. 7.678073e-04 9.593636e-04 8.608960e-04
## p.47.1. 6.571725e-04 8.091164e-04 7.309015e-04
## p.48.1. 1.731517e-04 2.741736e-04 2.209596e-04
## p.49.1. 5.705887e-04 7.723806e-04 6.680907e-04
## p.50.1. 6.340084e-04 9.435087e-04 7.811158e-04
## p.51.1. 1.250224e-04 2.347099e-04 1.759274e-04
## p.52.1. 9.767028e-05 1.986417e-04 1.435782e-04
## p.53.1. 2.362311e-04 3.254346e-04 2.794207e-04
## p.54.1. 1.663558e-04 2.688327e-04 2.143053e-04
## p.55.1. 1.970181e-04 3.039426e-04 2.473618e-04
## p.56.1. 4.027812e-04 7.881140e-04 5.795625e-04
## p.57.1. 6.524867e-05 1.308392e-04 9.527735e-05
## p.58.1. 9.776976e-05 1.966325e-04 1.427305e-04
## p.1.2. 5.728648e-04 9.569328e-04 7.527575e-04
## p.2.2. 3.123056e-03 5.480685e-03 4.216274e-03
## p.3.2. 2.541093e-03 5.475114e-03 3.866558e-03
## p.4.2. 4.208263e-03 9.728402e-03 6.682934e-03
## p.5.2. 2.895134e-04 4.590987e-04 3.693698e-04
## p.6.2. 7.331844e-04 1.608890e-03 1.132271e-03
## p.7.2. 4.092683e-04 6.048801e-04 5.023660e-04
## p.8.2. 2.526283e-04 3.372772e-04 2.932102e-04
## p.9.2. 5.295666e-04 6.648539e-04 5.950583e-04
## p.10.2. 2.685860e-04 3.966341e-04 3.292704e-04
## p.11.2. 5.237104e-04 7.242380e-04 6.199842e-04
## p.12.2. 7.005456e-04 9.082905e-04 8.012402e-04
## p.13.2. 1.686131e-04 2.596814e-04 2.117505e-04
## p.14.2. 2.635849e-04 3.571173e-04 3.084938e-04
## p.15.2. 3.337164e-04 4.713286e-04 3.995741e-04
## p.16.2. 2.867107e-04 4.029180e-04 3.419915e-04
## p.17.2. 4.445524e-04 5.513344e-04 4.963503e-04
## p.18.2. 2.296547e-04 3.016641e-04 2.644414e-04
## p.19.2. 6.254789e-04 7.694366e-04 6.955662e-04
## p.20.2. 7.360015e-04 8.813905e-04 8.075264e-04
## p.21.2. 1.070397e-03 1.331154e-03 1.196273e-03
## p.22.2. 3.378522e-04 4.574546e-04 3.951467e-04
## p.23.2. 5.307705e-04 6.458899e-04 5.867551e-04
## p.24.2. 1.182316e-03 1.436367e-03 1.305999e-03
## p.25.2. 1.887319e-04 2.622950e-04 2.238210e-04
## p.26.2. 2.510815e-04 3.579452e-04 3.020545e-04
## p.27.2. 6.735556e-04 8.139719e-04 7.421475e-04
## p.28.2. 3.431464e-04 4.244278e-04 3.826273e-04
## p.29.2. 1.833713e-04 2.698795e-04 2.246621e-04
## p.30.2. 7.133952e-05 1.561970e-04 1.097630e-04
## p.31.2. 2.369008e-04 3.130783e-04 2.737190e-04
## p.32.2. 9.349718e-05 2.439426e-04 1.599960e-04
## p.33.2. 1.686401e-04 2.420886e-04 2.035028e-04
## p.34.2. 9.115470e-04 1.135974e-03 1.020161e-03
## p.35.2. 7.112579e-04 8.920289e-04 7.991898e-04
## p.36.2. 1.328817e-04 2.090719e-04 1.691687e-04
## p.37.2. 3.138211e-04 4.002093e-04 3.557287e-04
## p.38.2. 2.923472e-04 3.938325e-04 3.413463e-04
## p.39.2. 2.516904e-04 3.414163e-04 2.947725e-04
## p.40.2. 4.555114e-04 7.013444e-04 5.714865e-04
## p.41.2. 3.858569e-04 5.151362e-04 4.478798e-04
## p.42.2. 9.426719e-04 1.148236e-03 1.042927e-03
## p.43.2. 1.684774e-04 2.623243e-04 2.128298e-04
## p.44.2. 6.175889e-05 1.436129e-04 9.848270e-05
## p.45.2. 3.025560e-04 3.898099e-04 3.447452e-04
## p.46.2. 7.650749e-04 9.560561e-04 8.578822e-04
## p.47.2. 6.547829e-04 8.062581e-04 7.282830e-04
## p.48.2. 1.724074e-04 2.730575e-04 2.200372e-04
## p.49.2. 5.684738e-04 7.696343e-04 6.656679e-04
## p.50.2. 6.316917e-04 9.402480e-04 7.783447e-04
## p.51.2. 1.244647e-04 2.337363e-04 1.751736e-04
## p.52.2. 9.722266e-05 1.978012e-04 1.429488e-04
## p.53.2. 2.352520e-04 3.241375e-04 2.782863e-04
## p.54.2. 1.656373e-04 2.677357e-04 2.134074e-04
## p.55.2. 1.961838e-04 3.027209e-04 2.463430e-04
## p.56.2. 4.012185e-04 7.853196e-04 5.774230e-04
## p.57.2. 6.493660e-05 1.302585e-04 9.484026e-05
## p.58.2. 9.732173e-05 1.957996e-04 1.421043e-04
## p.1.3. 5.707437e-04 9.536341e-04 7.500757e-04
## p.2.3. 3.114142e-03 5.466606e-03 4.204903e-03
## p.3.3. 2.533576e-03 5.461046e-03 3.855977e-03
## p.4.3. 4.196887e-03 9.706281e-03 6.666523e-03
## p.5.3. 2.883433e-04 4.573482e-04 3.679228e-04
## p.6.3. 7.305597e-04 1.603762e-03 1.128474e-03
## p.7.3. 4.076845e-04 6.026567e-04 5.004743e-04
## p.8.3. 2.515901e-04 3.359396e-04 2.920273e-04
## p.9.3. 5.275850e-04 6.624414e-04 5.928667e-04
## p.10.3. 2.674904e-04 3.950929e-04 3.279613e-04
## p.11.3. 5.217478e-04 7.216408e-04 6.177140e-04
## p.12.3. 6.980218e-04 9.051359e-04 7.984081e-04
## p.13.3. 1.678864e-04 2.586178e-04 2.108624e-04
## p.14.3. 2.625072e-04 3.557111e-04 3.072571e-04
## p.15.3. 3.323911e-04 4.695376e-04 3.980238e-04
## p.16.3. 2.855504e-04 4.013556e-04 3.406382e-04
## p.17.3. 4.428503e-04 5.492824e-04 4.944775e-04
## p.18.3. 2.287001e-04 3.004510e-04 2.633610e-04
## p.19.3. 6.231903e-04 7.667006e-04 6.930584e-04
## p.20.3. 7.333681e-04 8.783161e-04 8.046748e-04
## p.21.3. 1.066767e-03 1.326784e-03 1.192284e-03
## p.22.3. 3.365125e-04 4.557096e-04 3.936112e-04
## p.23.3. 5.287850e-04 6.435369e-04 5.845899e-04
## p.24.3. 1.178365e-03 1.431707e-03 1.301700e-03
## p.25.3. 1.879290e-04 2.612220e-04 2.228881e-04
## p.26.3. 2.500490e-04 3.565361e-04 3.008406e-04
## p.27.3. 6.711158e-04 8.111003e-04 7.394956e-04
## p.28.3. 3.417884e-04 4.227930e-04 3.811340e-04
## p.29.3. 1.825886e-04 2.687793e-04 2.237263e-04
## p.30.3. 7.100164e-05 1.555179e-04 1.092676e-04
## p.31.3. 2.359196e-04 3.118250e-04 2.726054e-04
## p.32.3. 9.306688e-05 2.429359e-04 1.593045e-04
## p.33.3. 1.679133e-04 2.410886e-04 2.026451e-04
## p.34.3. 9.083827e-04 1.132156e-03 1.016678e-03
## p.35.3. 7.087009e-04 8.889227e-04 7.963638e-04
## p.36.3. 1.322933e-04 2.081931e-04 1.684404e-04
## p.37.3. 3.125652e-04 3.986560e-04 3.543277e-04
## p.38.3. 2.911670e-04 3.923009e-04 3.399951e-04
## p.39.3. 2.506556e-04 3.400643e-04 2.935842e-04
## p.40.3. 4.537728e-04 6.988182e-04 5.693715e-04
## p.41.3. 3.843524e-04 5.132016e-04 4.461672e-04
## p.42.3. 9.394154e-04 1.144383e-03 1.039378e-03
## p.43.3. 1.677512e-04 2.612512e-04 2.119377e-04
## p.44.3. 6.146197e-05 1.429825e-04 9.803308e-05
## p.45.3. 3.013397e-04 3.882920e-04 3.433821e-04
## p.46.3. 7.623523e-04 9.527600e-04 8.548789e-04
## p.47.3. 6.524019e-04 8.034100e-04 7.256739e-04
## p.48.3. 1.716662e-04 2.719459e-04 2.191186e-04
## p.49.3. 5.663667e-04 7.668976e-04 6.632538e-04
## p.50.3. 6.293834e-04 9.369986e-04 7.755835e-04
## p.51.3. 1.239095e-04 2.327668e-04 1.744230e-04
## p.52.3. 9.677709e-05 1.969644e-04 1.423221e-04
## p.53.3. 2.342769e-04 3.228456e-04 2.771565e-04
## p.54.3. 1.649220e-04 2.666432e-04 2.125134e-04
## p.55.3. 1.953529e-04 3.015040e-04 2.453283e-04
## p.56.3. 3.996618e-04 7.825351e-04 5.752914e-04
## p.57.3. 6.462602e-05 1.296804e-04 9.440517e-05
## p.58.3. 9.687576e-05 1.949702e-04 1.414809e-04
## p.1.4. 5.686303e-04 9.503467e-04 7.474036e-04
## p.2.4. 3.105254e-03 5.452563e-03 4.193562e-03
## p.3.4. 2.526081e-03 5.447015e-03 3.845425e-03
## p.4.4. 4.185542e-03 9.684209e-03 6.650153e-03
## p.5.4. 2.871778e-04 4.556044e-04 3.664815e-04
## p.6.4. 7.279443e-04 1.598650e-03 1.124690e-03
## p.7.4. 4.061068e-04 6.004415e-04 4.985897e-04
## p.8.4. 2.505562e-04 3.346072e-04 2.908492e-04
## p.9.4. 5.256109e-04 6.600377e-04 5.906832e-04
## p.10.4. 2.663993e-04 3.935578e-04 3.266575e-04
## p.11.4. 5.197927e-04 7.190530e-04 6.154521e-04
## p.12.4. 6.955071e-04 9.019922e-04 7.955861e-04
## p.13.4. 1.671628e-04 2.575585e-04 2.099780e-04
## p.14.4. 2.614340e-04 3.543105e-04 3.060253e-04
## p.15.4. 3.310710e-04 4.677535e-04 3.964795e-04
## p.16.4. 2.843949e-04 3.997992e-04 3.392902e-04
## p.17.4. 4.411547e-04 5.472381e-04 4.926116e-04
## p.18.4. 2.277494e-04 2.992427e-04 2.622850e-04
## p.19.4. 6.209100e-04 7.639744e-04 6.905596e-04
## p.20.4. 7.307441e-04 8.752524e-04 8.018333e-04
## p.21.4. 1.063149e-03 1.322429e-03 1.188308e-03
## p.22.4. 3.351782e-04 4.539713e-04 3.920817e-04
## p.23.4. 5.268070e-04 6.411925e-04 5.824326e-04
## p.24.4. 1.174427e-03 1.427062e-03 1.297416e-03
## p.25.4. 1.871295e-04 2.601534e-04 2.219592e-04
## p.26.4. 2.490206e-04 3.551326e-04 2.996316e-04
## p.27.4. 6.686850e-04 8.082389e-04 7.368532e-04
## p.28.4. 3.404358e-04 4.211644e-04 3.796466e-04
## p.29.4. 1.818092e-04 2.676835e-04 2.227944e-04
## p.30.4. 7.066536e-05 1.548417e-04 1.087745e-04
## p.31.4. 2.349426e-04 3.105768e-04 2.714963e-04
## p.32.4. 9.263855e-05 2.419334e-04 1.586160e-04
## p.33.4. 1.671896e-04 2.400928e-04 2.017910e-04
## p.34.4. 9.052294e-04 1.128350e-03 1.013206e-03
## p.35.4. 7.061531e-04 8.858273e-04 7.935478e-04
## p.36.4. 1.317075e-04 2.073180e-04 1.677152e-04
## p.37.4. 3.113144e-04 3.971088e-04 3.529321e-04
## p.38.4. 2.899915e-04 3.907752e-04 3.386492e-04
## p.39.4. 2.496250e-04 3.387177e-04 2.924006e-04
## p.40.4. 4.520409e-04 6.963010e-04 5.672643e-04
## p.41.4. 3.828537e-04 5.112742e-04 4.444612e-04
## p.42.4. 9.361700e-04 1.140542e-03 1.035840e-03
## p.43.4. 1.670281e-04 2.601825e-04 2.110494e-04
## p.44.4. 6.116648e-05 1.423549e-04 9.758552e-05
## p.45.4. 3.001283e-04 3.867799e-04 3.420243e-04
## p.46.4. 7.596393e-04 9.494752e-04 8.518861e-04
## p.47.4. 6.500296e-04 8.005719e-04 7.230742e-04
## p.48.4. 1.709282e-04 2.708388e-04 2.182039e-04
## p.49.4. 5.642674e-04 7.641708e-04 6.608485e-04
## p.50.4. 6.270836e-04 9.337604e-04 7.728321e-04
## p.51.4. 1.233568e-04 2.318012e-04 1.736756e-04
## p.52.4. 9.633356e-05 1.961310e-04 1.416982e-04
## p.53.4. 2.333058e-04 3.215588e-04 2.760313e-04
## p.54.4. 1.642097e-04 2.655551e-04 2.116230e-04
## p.55.4. 1.945256e-04 3.002921e-04 2.443178e-04
## p.56.4. 3.981112e-04 7.797605e-04 5.731676e-04
## p.57.4. 6.431692e-05 1.291049e-04 9.397207e-05
## p.58.4. 9.643183e-05 1.941443e-04 1.408602e-04
## p.1.5. 5.665248e-04 9.470707e-04 7.447409e-04
## p.2.5. 3.096391e-03 5.438556e-03 4.182252e-03
## p.3.5. 2.518609e-03 5.433019e-03 3.834902e-03
## p.4.5. 4.174227e-03 9.662187e-03 6.633822e-03
## p.5.5. 2.860171e-04 4.538672e-04 3.650458e-04
## p.6.5. 7.253384e-04 1.593554e-03 1.120919e-03
## p.7.5. 4.045352e-04 5.982345e-04 4.967122e-04
## p.8.5. 2.495265e-04 3.332802e-04 2.896758e-04
## p.9.5. 5.236441e-04 6.576426e-04 5.885077e-04
## p.10.5. 2.653126e-04 3.920286e-04 3.253589e-04
## p.11.5. 5.178448e-04 7.164745e-04 6.131985e-04
## p.12.5. 6.930015e-04 8.988594e-04 7.927740e-04
## p.13.5. 1.664423e-04 2.565036e-04 2.090973e-04
## p.14.5. 2.603651e-04 3.529153e-04 3.047985e-04
## p.15.5. 3.297563e-04 4.659761e-04 3.949412e-04
## p.16.5. 2.832440e-04 3.982489e-04 3.379475e-04
## p.17.5. 4.394656e-04 5.452014e-04 4.907529e-04
## p.18.5. 2.268027e-04 2.980393e-04 2.612134e-04
## p.19.5. 6.186381e-04 7.612578e-04 6.880698e-04
## p.20.5. 7.281295e-04 8.721994e-04 7.990018e-04
## p.21.5. 1.059544e-03 1.318089e-03 1.184345e-03
## p.22.5. 3.338492e-04 4.522395e-04 3.905581e-04
## p.23.5. 5.248364e-04 6.388566e-04 5.802833e-04
## p.24.5. 1.170503e-03 1.422432e-03 1.293146e-03
## p.25.5. 1.863334e-04 2.590891e-04 2.210341e-04
## p.26.5. 2.479965e-04 3.537346e-04 2.984274e-04
## p.27.5. 6.662629e-04 8.053876e-04 7.342202e-04
## p.28.5. 3.390885e-04 4.195422e-04 3.781650e-04
## p.29.5. 1.810332e-04 2.665922e-04 2.218664e-04
## p.30.5. 7.033066e-05 1.541684e-04 1.082836e-04
## p.31.5. 2.339696e-04 3.093335e-04 2.703917e-04
## p.32.5. 9.221220e-05 2.409349e-04 1.579305e-04
## p.33.5. 1.664690e-04 2.391010e-04 2.009406e-04
## p.34.5. 9.020869e-04 1.124557e-03 1.009747e-03
## p.35.5. 7.036145e-04 8.827427e-04 7.907417e-04
## p.36.5. 1.311243e-04 2.064465e-04 1.669932e-04
## p.37.5. 3.100685e-04 3.955676e-04 3.515420e-04
## p.38.5. 2.888208e-04 3.892555e-04 3.373086e-04
## p.39.5. 2.485987e-04 3.373763e-04 2.912218e-04
## p.40.5. 4.503156e-04 6.937929e-04 5.651649e-04
## p.41.5. 3.813609e-04 5.093541e-04 4.427617e-04
## p.42.5. 9.329359e-04 1.136714e-03 1.032314e-03
## p.43.5. 1.663081e-04 2.591181e-04 2.101648e-04
## p.44.5. 6.087241e-05 1.417300e-04 9.714000e-05
## p.45.5. 2.989218e-04 3.852737e-04 3.406719e-04
## p.46.5. 7.569359e-04 9.462018e-04 8.489037e-04
## p.47.5. 6.476660e-04 7.977438e-04 7.204837e-04
## p.48.5. 1.701933e-04 2.697363e-04 2.172930e-04
## p.49.5. 5.621759e-04 7.614536e-04 6.584519e-04
## p.50.5. 6.247921e-04 9.305334e-04 7.700904e-04
## p.51.5. 1.228066e-04 2.308397e-04 1.729315e-04
## p.52.5. 9.589207e-05 1.953012e-04 1.410770e-04
## p.53.5. 2.323388e-04 3.202771e-04 2.749106e-04
## p.54.5. 1.635005e-04 2.644714e-04 2.107364e-04
## p.55.5. 1.937018e-04 2.990850e-04 2.433114e-04
## p.56.5. 3.965666e-04 7.769957e-04 5.710517e-04
## p.57.5. 6.400930e-05 1.285319e-04 9.354097e-05
## p.58.5. 9.598994e-05 1.933219e-04 1.402422e-04
## p.1.6. 5.644271e-04 9.438060e-04 7.420878e-04
## p.2.6. 3.087553e-03 5.424585e-03 4.170972e-03
## p.3.6. 2.511159e-03 5.419059e-03 3.824407e-03
## p.4.6. 4.162943e-03 9.640215e-03 6.617532e-03
## p.5.6. 2.848610e-04 4.521367e-04 3.636157e-04
## p.6.6. 7.227417e-04 1.588474e-03 1.117160e-03
## p.7.6. 4.029696e-04 5.960356e-04 4.948418e-04
## p.8.6. 2.485011e-04 3.319584e-04 2.885072e-04
## p.9.6. 5.216847e-04 6.552563e-04 5.863403e-04
## p.10.6. 2.642304e-04 3.905054e-04 3.240654e-04
## p.11.6. 5.159042e-04 7.139052e-04 6.109531e-04
## p.12.6. 6.905049e-04 8.957375e-04 7.899718e-04
## p.13.6. 1.657249e-04 2.554530e-04 2.082204e-04
## p.14.6. 2.593006e-04 3.515256e-04 3.035766e-04
## p.15.6. 3.284467e-04 4.642055e-04 3.934088e-04
## p.16.6. 2.820978e-04 3.967046e-04 3.366102e-04
## p.17.6. 4.377829e-04 5.431723e-04 4.889011e-04
## p.18.6. 2.258599e-04 2.968408e-04 2.601462e-04
## p.19.6. 6.163745e-04 7.585508e-04 6.855890e-04
## p.20.6. 7.255242e-04 8.691570e-04 7.961803e-04
## p.21.6. 1.055950e-03 1.313762e-03 1.180395e-03
## p.22.6. 3.325254e-04 4.505144e-04 3.890404e-04
## p.23.6. 5.228731e-04 6.365292e-04 5.781420e-04
## p.24.6. 1.166592e-03 1.417818e-03 1.288890e-03
## p.25.6. 1.855407e-04 2.580292e-04 2.201129e-04
## p.26.6. 2.469766e-04 3.523422e-04 2.972281e-04
## p.27.6. 6.638495e-04 8.025463e-04 7.315966e-04
## p.28.6. 3.377466e-04 4.179261e-04 3.766891e-04
## p.29.6. 1.802604e-04 2.655054e-04 2.209423e-04
## p.30.6. 6.999756e-05 1.534981e-04 1.077949e-04
## p.31.6. 2.330006e-04 3.080953e-04 2.692917e-04
## p.32.6. 9.178781e-05 2.399406e-04 1.572480e-04
## p.33.6. 1.657515e-04 2.381134e-04 2.000937e-04
## p.34.6. 8.989554e-04 1.120777e-03 1.006299e-03
## p.35.6. 7.010849e-04 8.796688e-04 7.879455e-04
## p.36.6. 1.305436e-04 2.055787e-04 1.662743e-04
## p.37.6. 3.088277e-04 3.940324e-04 3.501574e-04
## p.38.6. 2.876548e-04 3.877416e-04 3.359733e-04
## p.39.6. 2.475766e-04 3.360403e-04 2.900478e-04
## p.40.6. 4.485968e-04 6.912939e-04 5.630733e-04
## p.41.6. 3.798738e-04 5.074411e-04 4.410686e-04
## p.42.6. 9.297129e-04 1.132899e-03 1.028801e-03
## p.43.6. 1.655913e-04 2.580581e-04 2.092839e-04
## p.44.6. 6.057975e-05 1.411078e-04 9.669652e-05
## p.45.6. 2.977201e-04 3.837734e-04 3.393248e-04
## p.46.6. 7.542422e-04 9.429396e-04 8.459319e-04
## p.47.6. 6.453109e-04 7.949257e-04 7.179025e-04
## p.48.6. 1.694617e-04 2.686382e-04 2.163858e-04
## p.49.6. 5.600921e-04 7.587460e-04 6.560640e-04
## p.50.6. 6.225091e-04 9.273175e-04 7.673584e-04
## p.51.6. 1.222588e-04 2.298822e-04 1.721905e-04
## p.52.6. 9.545260e-05 1.944749e-04 1.404585e-04
## p.53.6. 2.313758e-04 3.190005e-04 2.737945e-04
## p.54.6. 1.627943e-04 2.633922e-04 2.098535e-04
## p.55.6. 1.928815e-04 2.978827e-04 2.423092e-04
## p.56.6. 3.950280e-04 7.742407e-04 5.689436e-04
## p.57.6. 6.370315e-05 1.279615e-04 9.311184e-05
## p.58.6. 9.555007e-05 1.925030e-04 1.396270e-04
## p.1.7. 5.623372e-04 9.405525e-04 7.394441e-04
## p.2.7. 3.078741e-03 5.410649e-03 4.159722e-03
## p.3.7. 2.503730e-03 5.405135e-03 3.813941e-03
## p.4.7. 4.151689e-03 9.618292e-03 6.601281e-03
## p.5.7. 2.837097e-04 4.504128e-04 3.621913e-04
## p.6.7. 7.201543e-04 1.583411e-03 1.113414e-03
## p.7.7. 4.014102e-04 5.938447e-04 4.929784e-04
## p.8.7. 2.474799e-04 3.306418e-04 2.873433e-04
## p.9.7. 5.197326e-04 6.528786e-04 5.841808e-04
## p.10.7. 2.631525e-04 3.889881e-04 3.227771e-04
## p.11.7. 5.139709e-04 7.113451e-04 6.087160e-04
## p.12.7. 6.880173e-04 8.926264e-04 7.871796e-04
## p.13.7. 1.650107e-04 2.544067e-04 2.073471e-04
## p.14.7. 2.582405e-04 3.501414e-04 3.023597e-04
## p.15.7. 3.271423e-04 4.624417e-04 3.918824e-04
## p.16.7. 2.809563e-04 3.951663e-04 3.352781e-04
## p.17.7. 4.361067e-04 5.411507e-04 4.870563e-04
## p.18.7. 2.249211e-04 2.956470e-04 2.590833e-04
## p.19.7. 6.141192e-04 7.558535e-04 6.831171e-04
## p.20.7. 7.229283e-04 8.661252e-04 7.933687e-04
## p.21.7. 1.052369e-03 1.309450e-03 1.176459e-03
## p.22.7. 3.312068e-04 4.487958e-04 3.875286e-04
## p.23.7. 5.209172e-04 6.342103e-04 5.760085e-04
## p.24.7. 1.162693e-03 1.413218e-03 1.284648e-03
## p.25.7. 1.847514e-04 2.569737e-04 2.191955e-04
## p.26.7. 2.459609e-04 3.509551e-04 2.960336e-04
## p.27.7. 6.614450e-04 7.997150e-04 7.289824e-04
## p.28.7. 3.364099e-04 4.163163e-04 3.752190e-04
## p.29.7. 1.794910e-04 2.644230e-04 2.200220e-04
## p.30.7. 6.966603e-05 1.528307e-04 1.073084e-04
## p.31.7. 2.320357e-04 3.068619e-04 2.681961e-04
## p.32.7. 9.136537e-05 2.389504e-04 1.565684e-04
## p.33.7. 1.650371e-04 2.371298e-04 1.992503e-04
## p.34.7. 8.958348e-04 1.117009e-03 1.002863e-03
## p.35.7. 6.985645e-04 8.766056e-04 7.851592e-04
## p.36.7. 1.299656e-04 2.047146e-04 1.655584e-04
## p.37.7. 3.075918e-04 3.925031e-04 3.487783e-04
## p.38.7. 2.864935e-04 3.862337e-04 3.346433e-04
## p.39.7. 2.465588e-04 3.347096e-04 2.888784e-04
## p.40.7. 4.468846e-04 6.888038e-04 5.609893e-04
## p.41.7. 3.783926e-04 5.055354e-04 4.393821e-04
## p.42.7. 9.265010e-04 1.129097e-03 1.025299e-03
## p.43.7. 1.648775e-04 2.570024e-04 2.084067e-04
## p.44.7. 6.028850e-05 1.404884e-04 9.625506e-05
## p.45.7. 2.965233e-04 3.822789e-04 3.379831e-04
## p.46.7. 7.515581e-04 9.396887e-04 8.429704e-04
## p.47.7. 6.429644e-04 7.921175e-04 7.153306e-04
## p.48.7. 1.687332e-04 2.675446e-04 2.154825e-04
## p.49.7. 5.580161e-04 7.560481e-04 6.536847e-04
## p.50.7. 6.202344e-04 9.241127e-04 7.646362e-04
## p.51.7. 1.217134e-04 2.289286e-04 1.714527e-04
## p.52.7. 9.501515e-05 1.936521e-04 1.398428e-04
## p.53.7. 2.304168e-04 3.177291e-04 2.726830e-04
## p.54.7. 1.620912e-04 2.623174e-04 2.089743e-04
## p.55.7. 1.920646e-04 2.966853e-04 2.413112e-04
## p.56.7. 3.934953e-04 7.714955e-04 5.668433e-04
## p.57.7. 6.339847e-05 1.273936e-04 9.268468e-05
## p.58.7. 9.511221e-05 1.916876e-04 1.390144e-04
## p.1.8. 5.602550e-04 9.373103e-04 7.368098e-04
## p.2.8. 3.069953e-03 5.396749e-03 4.148503e-03
## p.3.8. 2.496324e-03 5.391246e-03 3.803504e-03
## p.4.8. 4.140466e-03 9.596419e-03 6.585070e-03
## p.5.8. 2.825629e-04 4.486954e-04 3.607724e-04
## p.6.8. 7.175762e-04 1.578364e-03 1.109681e-03
## p.7.8. 3.998567e-04 5.916619e-04 4.911220e-04
## p.8.8. 2.464629e-04 3.293304e-04 2.861841e-04
## p.9.8. 5.177879e-04 6.505095e-04 5.820292e-04
## p.10.8. 2.620791e-04 3.874766e-04 3.214939e-04
## p.11.8. 5.120449e-04 7.087941e-04 6.064870e-04
## p.12.8. 6.855386e-04 8.895261e-04 7.843972e-04
## p.13.8. 1.642994e-04 2.533647e-04 2.064774e-04
## p.14.8. 2.571847e-04 3.487627e-04 3.011475e-04
## p.15.8. 3.258431e-04 4.606845e-04 3.903620e-04
## p.16.8. 2.798193e-04 3.936340e-04 3.339513e-04
## p.17.8. 4.344369e-04 5.391366e-04 4.852185e-04
## p.18.8. 2.239861e-04 2.944581e-04 2.580248e-04
## p.19.8. 6.118721e-04 7.531657e-04 6.806541e-04
## p.20.8. 7.203416e-04 8.631040e-04 7.905671e-04
## p.21.8. 1.048800e-03 1.305152e-03 1.172535e-03
## p.22.8. 3.298935e-04 4.470838e-04 3.860227e-04
## p.23.8. 5.189685e-04 6.318999e-04 5.738829e-04
## p.24.8. 1.158808e-03 1.408633e-03 1.280420e-03
## p.25.8. 1.839655e-04 2.559224e-04 2.182819e-04
## p.26.8. 2.449494e-04 3.495736e-04 2.948439e-04
## p.27.8. 6.590491e-04 7.968937e-04 7.263775e-04
## p.28.8. 3.350786e-04 4.147127e-04 3.737547e-04
## p.29.8. 1.787248e-04 2.633450e-04 2.191055e-04
## p.30.8. 6.933607e-05 1.521662e-04 1.068241e-04
## p.31.8. 2.310747e-04 3.056336e-04 2.671049e-04
## p.32.8. 9.094488e-05 2.379643e-04 1.558917e-04
## p.33.8. 1.643258e-04 2.361503e-04 1.984105e-04
## p.34.8. 8.927250e-04 1.113254e-03 9.994387e-04
## p.35.8. 6.960531e-04 8.735531e-04 7.823827e-04
## p.36.8. 1.293901e-04 2.038541e-04 1.648457e-04
## p.37.8. 3.063609e-04 3.909798e-04 3.474046e-04
## p.38.8. 2.853370e-04 3.847316e-04 3.333185e-04
## p.39.8. 2.455451e-04 3.333841e-04 2.877138e-04
## p.40.8. 4.451790e-04 6.863227e-04 5.589132e-04
## p.41.8. 3.769172e-04 5.036368e-04 4.377020e-04
## p.42.8. 9.233003e-04 1.125307e-03 1.021809e-03
## p.43.8. 1.641668e-04 2.559511e-04 2.075332e-04
## p.44.8. 5.999865e-05 1.398717e-04 9.581562e-05
## p.45.8. 2.953313e-04 3.807903e-04 3.366466e-04
## p.46.8. 7.488835e-04 9.364490e-04 8.400192e-04
## p.47.8. 6.406264e-04 7.893193e-04 7.127679e-04
## p.48.8. 1.680078e-04 2.664554e-04 2.145829e-04
## p.49.8. 5.559478e-04 7.533597e-04 6.513141e-04
## p.50.8. 6.179679e-04 9.209190e-04 7.619236e-04
## p.51.8. 1.211705e-04 2.279790e-04 1.707181e-04
## p.52.8. 9.457969e-05 1.928328e-04 1.392297e-04
## p.53.8. 2.294617e-04 3.164627e-04 2.715759e-04
## p.54.8. 1.613912e-04 2.612469e-04 2.080988e-04
## p.55.8. 1.912512e-04 2.954927e-04 2.403172e-04
## p.56.8. 3.919686e-04 7.687600e-04 5.647507e-04
## p.57.8. 6.309524e-05 1.268282e-04 9.225948e-05
## p.58.8. 9.467636e-05 1.908756e-04 1.384045e-04
The 95% confidence interval for \(\beta\) is (-0.004096, 0.0002897). Since \(\beta\) is not considered signifant in the 5% level, we can proceed with Model 2 instead, in which the overall time trend is removed from the model.
#Alter the data frame for visualizations
zip_data <- data.frame(zipcode=c(1:58), zip=rownames(Crime_count))
crime_viz_1 <- crime_data_1 %>%
colMeans() %>%
as.data.frame() %>%
rownames_to_column(var="zipcode") %>%
setNames(., c("zipcode", "delta")) %>%
filter(grepl("delta", zipcode)) %>%
mutate(zipcode = as.integer(gsub("\\D", "", zipcode))) %>%
left_join(zip_data, by="zipcode") %>%
mutate(id=as.character(zip))
zip <- readOGR(dsn="Boundaries - ZIP Codes", layer="geo_export_e1262361-5c82-45ee-8427-ef228e06dc4a")
## OGR data source with driver: ESRI Shapefile
## Source: "Boundaries - ZIP Codes", layer: "geo_export_e1262361-5c82-45ee-8427-ef228e06dc4a"
## with 61 features
## It has 4 fields
zip_df <- fortify(zip, region = "zip")
crime_model_1_delta <- left_join(zip_df, crime_viz_1, by="id")
head(crime_model_1_delta)
## long lat order hole piece id group zipcode delta
## 1 -87.62271 41.88884 1 FALSE 1 60601 60601.1 1 -0.003580096
## 2 -87.62232 41.88879 2 FALSE 1 60601 60601.1 1 -0.003580096
## 3 -87.62120 41.88866 3 FALSE 1 60601 60601.1 1 -0.003580096
## 4 -87.62064 41.88860 4 FALSE 1 60601 60601.1 1 -0.003580096
## 5 -87.62058 41.88859 5 FALSE 1 60601 60601.1 1 -0.003580096
## 6 -87.62055 41.88859 6 FALSE 1 60601 60601.1 1 -0.003580096
## zip
## 1 60601
## 2 60601
## 3 60601
## 4 60601
## 5 60601
## 6 60601
model_1_map <- ggplot()
model_1_map <- model_1_map + geom_polygon(data = crime_model_1_delta, aes(x=long, y=lat, group=group, fill=delta), color = "black", size=0.2)
model_1_map <- model_1_map + coord_map()
model_1_map <- model_1_map + labs(x="", y="", title=paste("Crime Differential Growth by Zipcode"), fill="total")
model_1_map <- model_1_map + cleanup
print(model_1_map)
i: zipcodes in Chicago (58 total in our data)
j: years (from 2010 to 2017)
y[i, j]: observed count of violent crimes in Chicago by year and zipcode
p[i, j]: probability of being a victim in location i and time j
n[i, j]: the population in zipcode i in time j.
time[j]: the time in years
\(\delta\)[i]: the spatial temporal interaction term
\(\sigma\)[i]: the income trend
income[i]: income for each zipcode in Chicago
This model resembles Model 2 but now our prediction of the crime rate accounts for the different income levels in various Chicago zipcodes.
\[y_{ij} \sim Bin(p_{ij}, n_{i})\] \[logit(p_{ij}) = \delta_{i} * time_{j} + \sigma_{i} * income_{i}\]
\[\delta_{i} \sim N(0, prec.delta)\] \[\sigma_{i} \sim N(0, prec.sigma)\]
\[prec.u \sim Gamma(0.5, 0.0005)\] \[prec.delta \sim Gamma(0.5, 0.0005)\] \[prec.sigma \sim Gamma(0.5, 0.0005)\]
#Specify the model
crime_model_3 <- "model{
#Data
for(i in 1:58) {
for(j in 1:8){
y[i,j] ~ dbin(p[i,j], n[i])
logit(p[i,j]) = delta[i] * time[j] + sigma[i] * income[i]
}
}
#Priors
for(i in 1:58){
delta[i] ~ dnorm(0, prec.delta)
sigma[i] ~ dnorm(0, prec.sigma)
}
#Hyperpriors
prec.delta ~ dgamma(0.5, 0.0005)
prec.sigma ~ dgamma(0.5, 0.0005)
}"
#set up an algorithm to simulate the posterior by
#combining the model (nba_model) and data (y)
#set the random number seed
crime_jags_3 <- jags.model(textConnection(crime_model_3), data=list(y=Crime_count, n=Population.vector, time=time, income=Model_Data$per_capita_income),
inits=list(.RNG.name="base::Wichmann-Hill", .RNG.seed=1989))
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 464
## Unobserved stochastic nodes: 118
## Total graph size: 2564
##
## Initializing model
#simulate a sample from the posterior
#note that we specify both mu and tau variables
crime_sim_3 <- coda.samples(crime_jags_3, variable.names=c("p", "delta", "sigma"), n.iter=40000)
#store the samples in a data frame:
crime_data_3 <- data.frame(crime_sim_3[[1]])
crime_viz_3_delta <- crime_data_3 %>%
colMeans() %>%
as.data.frame() %>%
rownames_to_column(var="zipcode") %>%
setNames(., c("zipcode", "delta")) %>%
filter(grepl("delta", zipcode)) %>%
mutate(zipcode = as.integer(gsub("\\D", "", zipcode))) %>%
left_join(zip_data, by="zipcode") %>%
mutate(id=as.character(zip))
crime_model_3_delta <- left_join(zip_df, crime_viz_3_delta, by="id")
head(crime_model_3_delta)
## long lat order hole piece id group zipcode delta
## 1 -87.62271 41.88884 1 FALSE 1 60601 60601.1 1 -0.1297566
## 2 -87.62232 41.88879 2 FALSE 1 60601 60601.1 1 -0.1297566
## 3 -87.62120 41.88866 3 FALSE 1 60601 60601.1 1 -0.1297566
## 4 -87.62064 41.88860 4 FALSE 1 60601 60601.1 1 -0.1297566
## 5 -87.62058 41.88859 5 FALSE 1 60601 60601.1 1 -0.1297566
## 6 -87.62055 41.88859 6 FALSE 1 60601 60601.1 1 -0.1297566
## zip
## 1 60601
## 2 60601
## 3 60601
## 4 60601
## 5 60601
## 6 60601
crime_viz_3_sigma <- crime_data_3 %>%
colMeans() %>%
as.data.frame() %>%
rownames_to_column(var="zipcode") %>%
setNames(., c("zipcode", "sigma")) %>%
filter(grepl("sigma", zipcode)) %>%
mutate(zipcode = as.integer(gsub("\\D", "", zipcode))) %>%
left_join(zip_data, by="zipcode") %>%
mutate(id=as.character(zip))
crime_model_3_sigma <- left_join(zip_df, crime_viz_3_sigma, by="id")
head(crime_model_3_sigma)
## long lat order hole piece id group zipcode sigma
## 1 -87.62271 41.88884 1 FALSE 1 60601 60601.1 1 0.002849706
## 2 -87.62232 41.88879 2 FALSE 1 60601 60601.1 1 0.002849706
## 3 -87.62120 41.88866 3 FALSE 1 60601 60601.1 1 0.002849706
## 4 -87.62064 41.88860 4 FALSE 1 60601 60601.1 1 0.002849706
## 5 -87.62058 41.88859 5 FALSE 1 60601 60601.1 1 0.002849706
## 6 -87.62055 41.88859 6 FALSE 1 60601 60601.1 1 0.002849706
## zip
## 1 60601
## 2 60601
## 3 60601
## 4 60601
## 5 60601
## 6 60601
model_3_map <- ggplot()
model_3_map <- model_3_map + geom_polygon(data = crime_model_3_delta, aes(x=long, y=lat, group=group, fill=delta), color = "black", size=0.2)
model_3_map <- model_3_map + coord_map()
model_3_map <- model_3_map + labs(x="", y="", title=paste("Crime Differential Growth by Zipcode"), fill="total")
model_3_map <- model_3_map + cleanup
print(model_3_map)
model_3_map <- ggplot()
model_3_map <- model_3_map + geom_polygon(data = crime_model_3_sigma, aes(x=long, y=lat, group=group, fill=sigma), color = "black", size=0.2)
model_3_map <- model_3_map + coord_map()
model_3_map <- model_3_map + labs(x="", y="", title=paste("Income Differential Growth by Zipcode"), fill="total")
model_3_map <- model_3_map + cleanup
print(model_3_map)
In locations where the income differential growth increases, the crime differential growth seems to decrease.
The running mean of \(\alpha\) appears to stabilize after 40,000 iterations.
running_mean_plot(crime_data_3$delta.1.)
running_mean_plot(crime_data_3$delta.3.)
running_mean_plot(crime_data_3$sigma.1.)
running_mean_plot(crime_data_3$sigma.3.)
i: zipcodes in Chicago (58 total in our data)
j: years (from 2010 to 2017)
y[i, j]: observed count of violent crimes in Chicago by year and zipcode
p[i, j]: probability of being a victim in location i and time j
n[i, j]: the population in zipcode i in time j.
\(\alpha\): the mean crime rate
u[i]: the unstructured spatial random effect for each zipcode in Chicago.
s[i]: the structured spatial random effect for each zipcode in Chicago.
time[j]: the time in years
\(\delta\)[i]: the spatial temporal interaction term
\(\sigma\)[i]: the income trend
income[i]: income for each zipcode in Chicago
e[i]: the weight in the relationship between crime rate and percentage of nonewhite people in Chicago
percent_nonwhite[i]: the percentage of nonwhite people in zipcode i.
This model resembles Model 3 but now our prediction of the crime rate accounts for both income levels and the percentage of nonwhites in the area.
\[y_{ij} \sim Bin(p_{ij}, n_{i})\] \[logit(p_{ij}) = \alpha + u_{i} + \delta_{i} * time_{j} + \sigma_{i} * income_{i} + \rho_{i} * percentNonwhite_{i}\]
\[\alpha \sim N(0, 1000^2)\] \[\beta \sim N(0, 1000^2)\] \[u_{i} \sim N(0, prec.u)\] \[\delta_{i} \sim N(0, prec.delta)\] \[\sigma_{i} \sim N(0, prec.sigma)\] \[\rho_{i} \sim N(0, 0.01)\]
\[prec.u \sim Gamma(0.5, 0.0005)\] \[prec.delta \sim Gamma(0.5, 0.0005)\] \[prec.sigma \sim Gamma(0.5, 0.0005)\]
#Specify the model
crime_model_4 <- "model{
#Data
for(i in 1:58) {
for(j in 1:8){
y[i,j] ~ dbin(p[i,j], n[i])
logit(p[i,j]) = delta[i] * time[j] + sigma[i]*income[i] + rho[i]*percent_nonwhite[i]
}
}
#Priors
for(i in 1:58){
delta[i] ~ dnorm(0, prec.delta)
sigma[i] ~ dnorm(0, prec.sigma)
rho[i] ~ dnorm(0, 0.01)
}
#Hyperpriors
prec.delta ~ dgamma(0.5, 0.0005)
prec.sigma ~ dgamma(0.5, 0.0005)
}"
#set up an algorithm to simulate the posterior by
#combining the model (nba_model) and data (y)
#set the random number seed
crime_jags_4 <- jags.model(textConnection(crime_model_4), data=list(y=Crime_count, n=Population.vector, time=time, income=Model_Data$per_capita_income, percent_nonwhite=100 - Model_Data$percent_white),
inits=list(.RNG.name="base::Wichmann-Hill", .RNG.seed=1989))
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 464
## Unobserved stochastic nodes: 176
## Total graph size: 3144
##
## Initializing model
#simulate a sample from the posterior
#note that we specify both mu and tau variables
crime_sim_4 <- coda.samples(crime_jags_4, variable.names=c("p", "delta", "sigma", "rho"), n.iter=40000)
#store the samples in a data frame:
crime_data_4 <- data.frame(crime_sim_4[[1]])
running_mean_plot(crime_data_4$delta.4.)
running_mean_plot(crime_data_4$sigma.4.)
running_mean_plot(crime_data_4$rho.4.)
crime_viz_4_delta <- crime_data_4 %>%
colMeans() %>%
as.data.frame() %>%
rownames_to_column(var="zipcode") %>%
setNames(., c("zipcode", "delta")) %>%
filter(grepl("delta", zipcode)) %>%
mutate(zipcode = as.integer(gsub("\\D", "", zipcode))) %>%
left_join(zip_data, by="zipcode") %>%
mutate(id=as.character(zip))
crime_model_4_delta <- left_join(zip_df, crime_viz_4_delta, by="id")
head(crime_model_4_delta)
## long lat order hole piece id group zipcode delta
## 1 -87.62271 41.88884 1 FALSE 1 60601 60601.1 1 -0.02008539
## 2 -87.62232 41.88879 2 FALSE 1 60601 60601.1 1 -0.02008539
## 3 -87.62120 41.88866 3 FALSE 1 60601 60601.1 1 -0.02008539
## 4 -87.62064 41.88860 4 FALSE 1 60601 60601.1 1 -0.02008539
## 5 -87.62058 41.88859 5 FALSE 1 60601 60601.1 1 -0.02008539
## 6 -87.62055 41.88859 6 FALSE 1 60601 60601.1 1 -0.02008539
## zip
## 1 60601
## 2 60601
## 3 60601
## 4 60601
## 5 60601
## 6 60601
crime_viz_4_sigma <- crime_data_4 %>%
colMeans() %>%
as.data.frame() %>%
rownames_to_column(var="zipcode") %>%
setNames(., c("zipcode", "sigma")) %>%
filter(grepl("sigma", zipcode)) %>%
mutate(zipcode = as.integer(gsub("\\D", "", zipcode))) %>%
left_join(zip_data, by="zipcode") %>%
mutate(id=as.character(zip))
crime_model_4_sigma <- left_join(zip_df, crime_viz_4_sigma, by="id")
head(crime_model_4_sigma)
## long lat order hole piece id group zipcode sigma
## 1 -87.62271 41.88884 1 FALSE 1 60601 60601.1 1 0.0001368507
## 2 -87.62232 41.88879 2 FALSE 1 60601 60601.1 1 0.0001368507
## 3 -87.62120 41.88866 3 FALSE 1 60601 60601.1 1 0.0001368507
## 4 -87.62064 41.88860 4 FALSE 1 60601 60601.1 1 0.0001368507
## 5 -87.62058 41.88859 5 FALSE 1 60601 60601.1 1 0.0001368507
## 6 -87.62055 41.88859 6 FALSE 1 60601 60601.1 1 0.0001368507
## zip
## 1 60601
## 2 60601
## 3 60601
## 4 60601
## 5 60601
## 6 60601
crime_viz_4_rho <- crime_data_4 %>%
colMeans() %>%
as.data.frame() %>%
rownames_to_column(var="zipcode") %>%
setNames(., c("zipcode", "rho")) %>%
filter(grepl("rho", zipcode)) %>%
mutate(zipcode = as.integer(gsub("\\D", "", zipcode))) %>%
left_join(zip_data, by="zipcode") %>%
mutate(id=as.character(zip))
crime_model_4_rho <- left_join(zip_df, crime_viz_4_rho, by="id")
head(crime_model_4_rho)
## long lat order hole piece id group zipcode rho
## 1 -87.62271 41.88884 1 FALSE 1 60601 60601.1 1 0.5684385
## 2 -87.62232 41.88879 2 FALSE 1 60601 60601.1 1 0.5684385
## 3 -87.62120 41.88866 3 FALSE 1 60601 60601.1 1 0.5684385
## 4 -87.62064 41.88860 4 FALSE 1 60601 60601.1 1 0.5684385
## 5 -87.62058 41.88859 5 FALSE 1 60601 60601.1 1 0.5684385
## 6 -87.62055 41.88859 6 FALSE 1 60601 60601.1 1 0.5684385
## zip
## 1 60601
## 2 60601
## 3 60601
## 4 60601
## 5 60601
## 6 60601
model_4_map <- ggplot()
model_4_map <- model_4_map + geom_polygon(data = crime_model_4_delta, aes(x=long, y=lat, group=group, fill=delta), color = "black", size=0.2)
model_4_map <- model_4_map + coord_map()
model_4_map <- model_4_map + labs(x="", y="", title=paste("Crime Differential Growth by Zipcode"), fill="total")
model_4_map <- model_4_map + cleanup
print(model_4_map)
model_4_map_sigma <- ggplot()
model_4_map_sigma <- model_4_map_sigma + geom_polygon(data = crime_model_4_sigma, aes(x=long, y=lat, group=group, fill=sigma), color = "black", size=0.2)
model_4_map_sigma <- model_4_map_sigma + coord_map()
model_4_map_sigma <- model_4_map_sigma + labs(x="", y="", title=paste("Income Differential Growth by Zipcode"), fill="total")
model_4_map_sigma <- model_4_map_sigma + cleanup
print(model_4_map_sigma)
model_4_map_rho <- ggplot()
model_4_map_rho <- model_4_map_rho + geom_polygon(data = crime_model_4_rho, aes(x=long, y=lat, group=group, fill=rho), color = "black", size=0.2)
model_4_map_rho <- model_4_map_rho + coord_map()
model_4_map_rho <- model_4_map_rho + labs(x="", y="", title=paste("Percentage Nonwhite Differential Growth by Zipcode"), fill="total")
model_4_map_rho <- model_4_map_rho + cleanup
print(model_4_map_rho)
#Number of zipcodes: 62; Number of races: 7
Police <- read.csv("Police_stops_2016_2017_sample.csv")
Police.r <- Police %>%
filter(!(RACE_CODE_CD=="I" | RACE_CODE_CD=="P" | RACE_CODE_CD=="WHI") ) %>%
mutate(RACE_CODE_CD = as.character(RACE_CODE_CD)) %>%
mutate(RACE_CODE_CD=replace(RACE_CODE_CD, RACE_CODE_CD=="WWH", "HISPANIC"))
vars = c("percent_white","percent_black","percent_asian","percent_hispanic")
Population.race <- Model_Data %>%
mutate_at(vars, function(x){trunc(Model_Data$total_population*x/100)}) %>%
select(zip, percent_asian,percent_black,percent_hispanic, percent_white) %>%
unique() %>%
remove_rownames %>%
column_to_rownames(var='zip')
head(Population.race)
## percent_asian percent_black percent_hispanic percent_white
## 60601 1751 1070 583 6131
## 60602 321 29 175 907
## 60603 264 35 44 538
## 60604 66 0 12 332
## 60605 4245 4494 1248 14483
## 60606 362 55 251 2091
valid_zips <- unique(rownames(Population.race))
Arrests <- Police.r %>%
filter(ZIP_CD %in% valid_zips) %>%
group_by(ZIP_CD, RACE_CODE_CD) %>%
dplyr::summarise(total.arrest = as.integer(sum(ENFORCEMENT_ACTION_TAKEN_I))) %>%
spread(key=ZIP_CD, value=total.arrest) %>%
as.data.frame() %>%
remove_rownames %>%
column_to_rownames(var='RACE_CODE_CD')
Arrests[is.na(Arrests)] <- as.integer(0)
head(Arrests)
## 60601 60602 60603 60604 60605 60606 60607 60608 60609 60610 60611
## API 0 0 0 0 0 0 0 0 1 1 0
## BLK 7 2 2 3 5 2 13 13 87 41 3
## HISPANIC 3 0 1 0 1 0 1 39 43 1 0
## WHT 2 0 0 0 4 0 5 7 8 3 0
## 60612 60613 60614 60615 60616 60617 60618 60619 60620 60621 60622
## API 0 1 0 0 3 0 0 0 1 0 0
## BLK 76 9 7 39 16 100 3 88 103 97 7
## HISPANIC 6 6 3 0 1 19 28 0 2 0 10
## WHT 9 7 4 0 3 3 11 0 3 0 6
## 60623 60624 60625 60626 60628 60629 60630 60631 60632 60633 60634
## API 0 2 1 0 4 0 0 0 0 0 0
## BLK 147 352 6 31 121 40 1 1 5 0 3
## HISPANIC 27 9 21 7 0 27 3 1 40 1 8
## WHT 7 12 4 2 1 4 7 3 5 0 4
## 60636 60637 60638 60639 60640 60641 60642 60643 60644 60645 60646
## API 1 0 0 0 1 0 0 0 0 0 0
## BLK 128 92 1 9 29 4 1 25 211 4 0
## HISPANIC 2 1 1 25 7 13 3 0 3 7 0
## WHT 2 0 4 4 7 8 1 2 9 3 1
## 60647 60649 60651 60652 60653 60654 60655 60656 60657 60659 60660
## API 0 0 0 0 0 0 0 0 0 3 2
## BLK 4 63 96 7 42 6 0 0 10 1 11
## HISPANIC 23 1 29 5 0 0 1 1 4 4 1
## WHT 3 1 6 2 0 1 2 0 2 6 2
## 60661 60707 60827
## API 0 0 0
## BLK 1 2 7
## HISPANIC 0 3 0
## WHT 1 2 0
Stops <- Police.r %>%
filter(ZIP_CD %in% valid_zips) %>%
group_by(ZIP_CD, RACE_CODE_CD) %>%
dplyr::summarise(total.stop = n()) %>%
spread(key=ZIP_CD, value=total.stop) %>%
as.data.frame() %>%
remove_rownames %>%
column_to_rownames(var='RACE_CODE_CD')
Stops[is.na(Stops)] <- 0
head(Stops)
## 60601 60602 60603 60604 60605 60606 60607 60608 60609 60610 60611
## API 0 0 0 0 0 0 0 0 2 2 0
## BLK 15 4 6 10 16 6 36 50 291 121 9
## HISPANIC 4 0 1 1 4 0 2 192 187 7 3
## WHT 4 0 2 1 5 0 11 27 36 11 4
## 60612 60613 60614 60615 60616 60617 60618 60619 60620 60621 60622
## API 1 1 1 0 4 1 0 1 2 0 1
## BLK 245 37 24 154 63 351 15 296 339 339 31
## HISPANIC 23 21 10 2 8 80 60 2 5 2 37
## WHT 24 24 12 4 8 13 22 1 5 2 21
## 60623 60624 60625 60626 60628 60629 60630 60631 60632 60633 60634
## API 1 4 6 2 4 3 8 0 1 0 0
## BLK 426 899 16 110 414 184 17 3 19 1 9
## HISPANIC 134 43 50 44 8 130 29 2 165 2 32
## WHT 21 71 32 27 9 29 64 9 28 0 23
## 60636 60637 60638 60639 60640 60641 60642 60643 60644 60645 60646
## API 3 2 2 0 5 0 0 0 1 3 1
## BLK 389 299 9 37 111 13 13 86 626 38 1
## HISPANIC 12 1 12 107 23 51 10 1 23 17 3
## WHT 6 5 18 17 25 25 5 9 25 13 4
## 60647 60649 60651 60652 60653 60654 60655 60656 60657 60659 60660
## API 1 4 0 0 1 0 0 1 0 8 5
## BLK 19 237 320 32 203 18 2 0 46 24 32
## HISPANIC 87 3 95 11 2 0 4 1 16 31 14
## WHT 8 5 41 10 4 7 7 5 19 23 7
## 60661 60707 60827
## API 1 0 0
## BLK 3 3 31
## HISPANIC 0 7 1
## WHT 4 2 0
total_pop <- Model_Data %>% filter(Year==2010) %>% select(total_population)
population_by_ethnic <- Model_Data %>%
filter(Year==2010) %>%
mutate_at(c("percent_asian", "percent_black", "percent_hispanic", "percent_white"), function(x){x* total_pop$total_population/100}) %>%
mutate(zip = as.numeric(zip)) %>%
select(percent_asian, percent_black, percent_hispanic, percent_white, zip)
population_by_ethnic
## percent_asian percent_black percent_hispanic percent_white zip
## 1 1751.76 1070.52 583.92 6131.16 60601
## 2 321.86 29.26 175.56 907.06 60602
## 3 264.60 35.28 44.10 538.02 60603
## 4 66.40 0.00 12.45 332.00 60604
## 5 4245.24 4494.96 1248.60 14483.76 60605
## 6 362.57 55.78 251.01 2091.75 60606
## 7 4775.46 4021.44 2010.72 13321.02 60607
## 8 7879.40 12607.04 45700.52 11819.10 60608
## 9 2510.92 17576.44 30758.77 10671.41 60609
## 10 2711.80 6585.80 2324.40 26343.20 60610
## 11 4458.30 1188.88 1188.88 22291.50 60611
## 12 1442.08 21991.72 5047.28 7210.40 60612
## 13 2972.16 2972.16 4953.60 37647.36 60613
## 14 3418.15 2734.52 4785.41 56057.66 60614
## 15 2852.15 24447.00 2037.25 10186.25 60615
## 16 19301.49 11877.84 5444.01 11877.84 60616
## 17 0.00 45054.90 30870.95 5840.45 60617
## 18 3896.52 1948.26 45784.11 42861.72 60618
## 19 0.00 62974.34 649.22 649.22 60619
## 20 0.00 71055.41 732.53 732.53 60620
## 21 0.00 32263.68 336.08 336.08 60621
## 22 1628.16 4341.76 15196.16 32020.48 60622
## 23 0.00 28953.72 53649.54 1703.16 60623
## 24 0.00 37318.85 785.66 785.66 60624
## 25 11055.24 3158.64 30007.08 33165.72 60625
## 26 3574.27 12765.25 11744.03 21445.62 60626
## 27 0.00 67407.40 2151.30 1434.20 60628
## 28 1138.33 26181.59 76268.11 10244.97 60629
## 29 6367.46 578.86 15629.22 34152.74 60630
## 30 577.98 288.99 3467.88 24275.16 60631
## 31 4542.25 1816.90 75401.35 9084.50 60632
## 32 134.74 3368.50 5659.08 4311.68 60633
## 33 2908.72 727.18 25451.30 42903.62 60634
## 34 0.00 37522.92 1197.54 399.18 60636
## 35 1930.80 37167.90 965.40 7723.20 60637
## 36 560.33 1680.99 22973.53 30257.82 60638
## 37 0.00 13790.10 70789.18 6435.38 60639
## 38 8367.45 11585.70 8367.45 34113.45 60640
## 39 2829.20 1414.60 38194.20 26877.40 60641
## 40 915.60 2014.32 3845.52 10804.08 60642
## 41 0.00 37321.92 1555.08 11922.28 60643
## 42 0.00 44773.14 1428.93 952.62 60644
## 43 7036.95 7506.08 10320.86 20641.72 60645
## 44 2831.40 283.14 3397.68 21518.64 60646
## 45 2614.38 5228.76 45315.92 33115.48 60647
## 46 0.00 42717.36 908.88 908.88 60649
## 47 0.00 38524.50 20179.50 1834.50 60651
## 48 0.00 20408.16 14880.95 6802.72 60652
## 49 310.86 28288.26 621.72 1243.44 60653
## 50 1588.50 794.25 1429.65 11596.05 60654
## 51 281.73 1972.11 1690.38 24228.78 60655
## 52 2197.44 274.68 3296.16 21150.36 60656
## 53 4766.16 2042.64 4766.16 55151.28 60657
## 54 12690.24 1982.85 7931.40 15862.80 60659
## 55 5467.80 6309.00 7570.80 21450.60 60660
## 56 1635.92 446.16 446.16 4759.04 60661
## 57 1281.36 2989.84 12386.48 24772.96 60707
## 58 0.00 26263.24 1141.88 856.41 60827
print(colSums(population_by_ethnic))
## percent_asian percent_black percent_hispanic percent_white
## 152464.1 885195.5 785951.5 893236.1
## zip
## 3516807.0
data_stops = data.frame(race=c("Asian Pacific Islander", "Black", "Hispanic", "White"), stop_rate=(rowSums(Stops)/colSums(population_by_ethnic %>% select(-zip))))
rownames(data_stops) <- c()
head(data_stops %>% arrange(desc(stop_rate)))
## race stop_rate
## 1 Black 0.0080750524
## 2 Hispanic 0.0023182092
## 3 White 0.0009448790
## 4 Asian Pacific Islander 0.0005443903
Blacks and Hispanics are stopped at higher rate than Whites and Asians given their population.
l <- population_by_ethnic %>%
mutate(BLK=percent_black, HISPANIC=percent_hispanic, API=percent_asian, WHT=percent_white) %>%
select(BLK, HISPANIC, API, WHT, zip) %>%
gather(key=RACE_CODE_CD, value=population, - zip)
data <- Police.r %>%
filter(ZIP_CD %in% valid_zips) %>%
group_by(ZIP_CD, RACE_CODE_CD) %>%
dplyr::summarise(total.stop = n()) %>%
inner_join(l, by=c("ZIP_CD" = 'zip', 'RACE_CODE_CD')) %>%
mutate(race=RACE_CODE_CD)
ggplot(data, aes(x=population, y=total.stop, colour=RACE_CODE_CD)) + geom_point() + geom_smooth() +
labs(x="Population Count", y="Number of Traffic Stops", title="Traffic Stops By Population", subtitle="Visualization of number of traffic stops as the population of each ethnic group increases for each zipcode") +
guides(colour=guide_legend("Race")) +
theme(plot.title = element_text(hjust=0.5))
## `geom_smooth()` using method = 'loess'
ggplot(data, aes(x=population, y=total.stop, colour=RACE_CODE_CD)) + geom_point() + geom_smooth() +
labs(x="Percent of Population", y="Number of Traffic Stops", fill="Race", title="Traffic Stops By Percentage of Population") +
guides(colour=guide_legend("Race")) +
theme(plot.title = element_text(hjust=0.5))
## `geom_smooth()` using method = 'loess'
Data is only from January 2016 - February 2017 (13/12)
i: ethnic group
j: zipcode locations in Chicago
y[i, j]: observed number of arrests of ethnic group i in area j
n[i, j]: the number of traffic stops in area j for ethnic group i
\(\mu\): grand mean of arrests
\(\alpha\)[i]: the random effect for each ethnic group i
\[y_{ij} \sim Pois(\frac{13}{12} n_{ij} \theta_{ij})\] \[log(\theta_{ij}) = \mu + \alpha_{i} \]
\[\mu \sim N(0, 0.000001)\] \[\alpha_{i} \sim N(0, 0.000001)\]
log(0.28603)
## [1] -1.251659
log(0.28603+0.05911) - log(0.28603)
## [1] 0.1878534
log(0.28603-0.05712) - log(0.28603)
## [1] -0.2227678
log(0.28603-0.12381) - log(0.28603)
## [1] -0.5671433
police_model_0 <- "model{
#Data
for (i in 1:4) {
for(j in 1:58){
y[i,j] ~ dpois((13/12) * n[i,j] * theta[i,j])
log(theta[i,j]) <- mu + alpha[i]
}
}
#Priors
mu ~ dnorm(0, 1.0E-06)
for (i in 1:4){
alpha[i] ~ dnorm(0, 0.005)
}
}"
#initial value set
inits = list(mu=-1.251659, alpha=c(0, 0.1878534, -0.2227678, -0.5671433))
police_jags_0 <- jags.model(textConnection(police_model_0), data=list(y=Arrests, n=Stops),
inits=inits)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 232
## Unobserved stochastic nodes: 5
## Total graph size: 597
##
## Initializing model
#simulate a sample from the posterior
#note that we specify both mu and tau variables
police_sim_0 <- coda.samples(police_jags_0, variable.names=c("mu", "alpha"), n.iter=40000)
#store the samples in a data frame:
police_data_0 <- data.frame(police_sim_0[[1]])
#estimate initial values
data <- Police.r %>%
filter(ZIP_CD %in% valid_zips) %>%
group_by(ZIP_CD, RACE_CODE_CD) %>%
dplyr::summarise(total.arrest = as.integer(sum(ENFORCEMENT_ACTION_TAKEN_I)), total.stop = n())
head(data)
## # A tibble: 6 x 4
## # Groups: ZIP_CD [3]
## ZIP_CD RACE_CODE_CD total.arrest total.stop
## <int> <chr> <int> <int>
## 1 60601 BLK 7 15
## 2 60601 HISPANIC 3 4
## 3 60601 WHT 2 4
## 4 60602 BLK 2 4
## 5 60603 BLK 2 6
## 6 60603 HISPANIC 1 1
freq_model <- lm(total.arrest ~ RACE_CODE_CD * total.stop, data = data)
summary(freq_model)
##
## Call:
## lm(formula = total.arrest ~ RACE_CODE_CD * total.stop, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23.097 -1.148 -0.147 1.273 46.686
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.08564 1.60786 -0.053 0.95758
## RACE_CODE_CDBLK -4.88013 1.84214 -2.649 0.00875 **
## RACE_CODE_CDHISPANIC 0.54742 1.85234 0.296 0.76791
## RACE_CODE_CDWHT 1.12422 1.95807 0.574 0.56655
## total.stop 0.28603 0.49111 0.582 0.56098
## RACE_CODE_CDBLK:total.stop 0.05911 0.49113 0.120 0.90433
## RACE_CODE_CDHISPANIC:total.stop -0.05712 0.49135 -0.116 0.90757
## RACE_CODE_CDWHT:total.stop -0.12381 0.49393 -0.251 0.80235
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.55 on 189 degrees of freedom
## Multiple R-squared: 0.9789, Adjusted R-squared: 0.9781
## F-statistic: 1250 on 7 and 189 DF, p-value: < 2.2e-16
Data is only from January 2016 - February 2017 (13/12)
i: ethnic group
j: zipcode locations in Chicago
y[i, j]: observed number of arrests of ethnic group i in area j
n[i, j]: the number of traffic stops in area j for ethnic group i
\(\mu\): grand mean of arrests
\(\alpha\)[i]: the random effect for each ethnic group i
\(\beta\)[j]: the random effect for each zipcode location j
\(\epsilon\)[i, j]: the random error for each ethnic group i and location j
\[y_{ij} \sim Pois(\frac{13}{12} n_{ij} \theta_{ij})\] \[log(\theta_{ij}) = \mu + \alpha_{i} + \beta_{j} + \epsilon_{ij}\]
\[\beta_{0} \sim N(0, 0.000001)\] \[\alpha_{i} \sim N(0, 0.000001)\] \[\beta_{j} \sim N(0, b^2)\] \[\epsilon_{ij} \sim N(0, e^2)\]
\[b \sim Gamma(0.5, 0.0005)\] \[e \sim Gamma(0.5, 0.0005)\]
police_model_1 <- "model{
#Data
for (i in 1:4) {
for(j in 1:58){
y[i,j] ~ dpois((13/12)*n[i,j]*theta[i,j])
log(theta[i,j]) <- mu + alpha[i] + beta[j] + epsilon[i,j]
}
}
#Priors
mu ~ dnorm(0,1.0E-06)
for (i in 1:4){
alpha[i] ~ dnorm(0, 0.005)
}
for (j in 1:58){
beta[j] ~ dnorm(0, 1/b^2)
}
for (i in 1:4){
for (j in 1:58){
epsilon[i,j] ~ dnorm(0, 1/e^2)
}
}
b ~ dgamma(0.5, 0.0005)
e ~ dgamma(0.5, 0.0005)
}"
#initial value set
inits = list(mu=-1.251659, alpha=c(0, 0.1878534, -0.2227678, -0.5671433), tau_beta=170, tau_epsilon=170, beta=rep(0.03, 58), epsilon=matrix(rep(0.01, 58*4), ncol=58))
police_jags_1 <- jags.model(textConnection(police_model_1), data=list(y=Arrests, n=Stops),
inits=inits)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 232
## Unobserved stochastic nodes: 297
## Total graph size: 1471
##
## Initializing model
#simulate a sample from the posterior
#note that we specify both mu and tau variables
police_sim_1 <- coda.samples(police_jags_1, variable.names=c("mu", "alpha", "beta", "sd_beta", "sd_epsilon", "epsilon"), n.iter=40000)
#store the samples in a data frame:
police_data_1 <- data.frame(police_sim_1[[1]])
running_mean_plot(police_data_1$beta.1.)
running_mean_plot(police_data_1$mu)
running_mean_plot(police_data_1$alpha.2.)
police_viz <- police_data_1 %>%
colMeans() %>%
as.data.frame() %>%
rownames_to_column(var="zipcode") %>%
setNames(., c("zipcode", "beta")) %>%
filter(grepl("beta", zipcode)) %>%
mutate(zipcode = as.integer(gsub("\\D", "", zipcode))) %>%
mutate(mu = mean(police_data_1$mu)) %>%
mutate(rate = exp(mu + beta)) %>%
left_join(zip_data, by="zipcode") %>%
mutate(id=as.character(zip))
police_model_beta <- left_join(zip_df, police_viz, by="id")
head(police_model_beta)
## long lat order hole piece id group zipcode beta
## 1 -87.62271 41.88884 1 FALSE 1 60601 60601.1 1 0.05094085
## 2 -87.62232 41.88879 2 FALSE 1 60601 60601.1 1 0.05094085
## 3 -87.62120 41.88866 3 FALSE 1 60601 60601.1 1 0.05094085
## 4 -87.62064 41.88860 4 FALSE 1 60601 60601.1 1 0.05094085
## 5 -87.62058 41.88859 5 FALSE 1 60601 60601.1 1 0.05094085
## 6 -87.62055 41.88859 6 FALSE 1 60601 60601.1 1 0.05094085
## mu rate zip
## 1 -0.2179287 0.8462099 60601
## 2 -0.2179287 0.8462099 60601
## 3 -0.2179287 0.8462099 60601
## 4 -0.2179287 0.8462099 60601
## 5 -0.2179287 0.8462099 60601
## 6 -0.2179287 0.8462099 60601
model_1_map_beta <- ggplot()
model_1_map_beta <- model_1_map_beta + geom_polygon(data = police_model_beta, aes(x=long, y=lat, group=group, fill=rate), color = "black", size=0.2)
model_1_map_beta <- model_1_map_beta + coord_map()
model_1_map_beta <- model_1_map_beta + labs(x="", y="", title=paste("Rate at Which Traffic Stops Result in Arrests by Zipcode"), fill="total")
model_1_map_beta <- model_1_map_beta + cleanup
print(model_1_map_beta)
police_viz1 <- police_data_1 %>%
colMeans() %>%
as.data.frame() %>%
rownames_to_column(var="zipcode") %>%
setNames(., c("zipcode", "alpha")) %>%
filter(grepl("alpha", zipcode)) %>%
mutate(mu = mean(police_data_1$mu)) %>%
mutate(rate = exp(mu + alpha))
head(police_viz1)
## zipcode alpha mu rate
## 1 alpha.1. -1.255402 -0.2179287 0.2291610
## 2 alpha.2. -1.097255 -0.2179287 0.2684249
## 3 alpha.3. -1.279364 -0.2179287 0.2237350
## 4 alpha.4. -1.344437 -0.2179287 0.2096395